SolarRadiation.f90 Source File
Manage solar radiation with arbitrary time cumulation
Source Code
!! Manage solar radiation with arbitrary time cumulation !|author: <a href="mailto:giovanni.ravazzani@polimi.it">Giovanni Ravazzani</a> ! license: <a href="http://www.gnu.org/licenses/">GPL</a> ! !### History ! ! current version 1.1 - 16th June 2021 ! ! | version | date | comment | ! |----------|-------------|----------| ! | 1.0 | 19/May/2021 | Original code | ! | 1.1 | 26/Jun/2021 | added net radiation computation | ! !### License ! license: GNU GPL <http://www.gnu.org/licenses/> ! !### Module Description ! Module to manage solar radiation with arbitrary time cumulation ! MODULE SolarRadiation ! Modules used: USE DataTypeSizes, ONLY : & ! Imported Parameters: float, short USE Loglib, ONLY : & !Imported routines: Catch USE ObservationalNetworks, ONLY : & !Imported types: ObservationalNetwork, & !Imported routines: ReadMetadata, ReadData, & CopyNetwork, & !Imported operators: ASSIGNMENT (=) USE Chronos, ONLY : & !Imported types: DateTime, & !Imported variables: timeString, & !imported routines: DayOfYear, GetHour, & !Imported operations: OPERATOR (<), & OPERATOR (>), & OPERATOR (>=), & OPERATOR (<=), & OPERATOR (+), & OPERATOR (==), & ASSIGNMENT (=) USE GridLib, ONLY : & !Imported types: grid_real, grid_integer, & !Imported routines: NewGrid, GetDtGrid, & GridDestroy, ExportGrid, & SetLongName, SetStandardName, & SetUnits, SetCurrentTime, & !Imported parameters: NET_CDF, ESRI_ASCII, & ESRI_BINARY USE GridOperations, ONLY : & !Imported routines: GridByIni, GridResample, & CRSisEqual, GridConvert, & GetXY, GetIJ, & !Imported operations: ASSIGNMENT (=) USE Units, ONLY : & !imported parameters: degToRad, pi, hour, stefanBoltzman USE IniLib, ONLY : & !Imported types: IniList, & !Imported routines: IniReadInt, IniReadReal, & KeyIsPresent, IniReadString USE Utilities, ONLY: & !Imported routines: GetUnit USE MeteoUtilities, ONLY: & !Imported routines: ReadField, Interpolate, & Offset, Scale USE FileSys, ONLY: & !Imported routines: DirExists, FileExists, & FileDelete, CurrentDir, & GetOS, DirNew, & !Imported parameters: WIN32 USE Morphology, ONLY: & !imported routines: DeriveAspect, DeriveSlope USE GeoLib, ONLY: & !imported type Coordinate, & !Imported variables: point1, point2, & !Imported routines: DecodeEPSG, Convert, & Distance, & !Imported types: CRS, & !Imported operations: OPERATOR (==) USE DomainProperties, ONLY : & !imported variables latCentroid USE StringManipulation, ONLY: & !Imported routines: ReplaceChar, ToString IMPLICIT NONE ! Global (i.e. public) Declarations: INTEGER (KIND = short) :: dtRadiation = 0!!cumulation deltat TYPE (ObservationalNetwork) :: radiometers !!radiation stations network TYPE (grid_real) :: radiation !![W/m2] TYPE (grid_real) :: netRadiation !!incoming and outcoming shortwave and longwave radiation [W/m2] TYPE (grid_real) :: netRadiationFAO !!net radiation of FAO reference grass with albedo = 0.23 !Private declarations !define type for viewing angle: the maximum angle of terrain obstruction TYPE ViewingAngle REAL (KIND = float) :: angle (16) END TYPE ViewingAngle TYPE(DateTime), PRIVATE :: timeNew !!time when new data must be read INTEGER (KIND = short), PRIVATE :: fileunit !! unit of file containing data INTEGER (KIND = short), PRIVATE :: interpolationMethod !!method to spatial interpolate site data INTEGER (KIND = short), PRIVATE :: interpolationMethod_assignment !!method to assign spatial interpolation ! !!1 = one method for the entire domain, 2 = a map with interpolation method codes TYPE (grid_integer), PRIVATE :: interpolationMethod_map INTEGER (KIND = short), PRIVATE :: interpolationMethod_vector (3) !!defines active interpolation methods TYPE (grid_real), PRIVATE :: interpolatedMap (3) !!1 = map for thiessen, 2 = map for idw, 3 = map for kriging REAL (KIND = float), PRIVATE :: cellsizeInterpolation !!spatial resolution of interpolated grid INTEGER (KIND = short), PRIVATE :: neighbors !!number of closest data to use for interpolation TYPE(grid_real), PRIVATE :: grid_devst !!standard deviation of kriging interpolation CHARACTER (LEN = 300), PRIVATE :: fileGrid !!file containing grids INTEGER (KIND = short), PRIVATE :: dtGrid !!dt of imported field INTEGER (KIND = short), PRIVATE :: elevationDrift !! 1 = use elevation to modify interpolated data INTEGER (KIND = short), PRIVATE :: timezone !!local timezone TYPE(grid_real), PRIVATE :: shadowGrid !!shadow grid 1 = shadow, 0 = no shadow TYPE (ViewingAngle), ALLOCATABLE, PRIVATE :: viewangle (:,:) !!raster of maximum viewing angles REAL (KIND = float), PRIVATE :: azimuth (16) TYPE(grid_real), PRIVATE :: slope !!terreain slope to be used to modify interpolated data [radians] TYPE(grid_real), PRIVATE :: aspect !!terreain aspect to be used to modify interpolated data [radians] INTEGER (KIND = short), PRIVATE :: i, j !!loop index REAL (KIND = float), PRIVATE, PARAMETER :: MISSING_DEF_REAL = -9999.9 REAL (KIND = float), PRIVATE :: scale_factor !!scale factor to apply to final map REAL (KIND = float), PRIVATE :: offset_value !!offset to apply to final map REAL (KIND = float), PRIVATE :: valid_prcn !!when data from several time steps are read, !!this is the minimum percentage (0-1) of valid !!data that must be prresent to consider !!valid the aggregated value. !used by inverse distance weighting interpolation REAL (KIND = float), PRIVATE :: idw_power !! power to be used with IDW !used by kriging INTEGER (KIND = short), PRIVATE :: krige_var !! when set to 1 a grid of kriging variance is generated and exported if export option is active, default to 0 INTEGER (KIND = short), PRIVATE :: krige_anisotropy !! considers anisotropy, default = 0 excludes anisotropy INTEGER (KIND = short), PRIVATE :: krige_varmodel !! 1 = spherical, 2 = exponential, 3 = gaussian, 0 = automatic fitting. default to 2 INTEGER (KIND = short), PRIVATE :: krige_lags !! number of lags for semivariogram. if undefined or set to 0 default to 15 REAL (KIND = float) , PRIVATE :: krige_maxlag !! maximum distance to be considered for semivariogram assessment. If undefined or set to 0, it is computed automatically INTEGER (KIND = short), PRIVATE :: export !!activates grid exporting CHARACTER (LEN = 1000), PRIVATE :: export_path !!folder where to put exported grids CHARACTER (LEN = 1000), PRIVATE :: export_file !!name of exported file CHARACTER (LEN = 1000), PRIVATE :: export_file_net !!name of exported net radiation file CHARACTER (LEN = 1000), PRIVATE :: export_file_var !!name of exported variance grid TYPE (DateTime), PRIVATE :: export_start !! time and date to start exporting TYPE (DateTime), PRIVATE :: export_stop !! time and date to stop exporting INTEGER (KIND = short), PRIVATE :: export_dt !! time between two exportations TYPE (DateTime), PRIVATE :: timeNewExport !! time when new exporting must occur INTEGER (KIND = short), PRIVATE :: export_format !! 1 = esri_ascii, 2 = esri_binary, 3 = netcdf INTEGER (KIND = short), PRIVATE :: export_epsg !! coordinate reference system of exported grid TYPE (grid_real), PRIVATE :: exportedGrid, exportedGridVar, exportedGridNet TYPE (CRS), PRIVATE :: exportGridMapping LOGICAL, PRIVATE :: needConversion !Public routines PUBLIC :: SolarRadiationInit PUBLIC :: SolarRadiationRead !Local routines PRIVATE :: SetSpecificProperties PRIVATE :: SkyView PRIVATE :: CastShadow PRIVATE :: SolarDeclination PRIVATE :: SunElevationAngle PRIVATE :: SolarHourAngle PRIVATE :: SolarAzimuth PRIVATE :: OpticalDepth PRIVATE :: ComputeNetRadiation !======= CONTAINS !======= ! Define procedures contained in this module. !============================================================================== !| Description: ! Initialize solar radiation SUBROUTINE SolarRadiationInit & ! ( ini, mask, dtMeteo, tstart, dem, dem_loaded, albedo_loaded, & dtTemperature, dtRelHumidity ) IMPLICIT NONE TYPE (IniList), INTENT(IN) :: ini TYPE (grid_integer), INTENT(IN) :: mask !!defines interpolation extent INTEGER (KIND = short), INTENT(IN) :: dtMeteo !! deltat of meteo data reading TYPE (DateTime), INTENT(IN) :: tstart !!initial time TYPE(grid_real), INTENT(in) :: dem !!digital elevation model to be used to modify interpolated data LOGICAL , INTENT (in) :: dem_loaded !! true if dem has been loaded LOGICAL , INTENT (in) :: albedo_loaded !! true if dem has been loaded INTEGER (KIND = short), INTENT(IN) :: dtTemperature !! delta time temperature INTEGER (KIND = short), INTENT(IN) :: dtRelHumidity !! delta time relative humidity !local declarations CHARACTER (LEN = 1000) :: filename TYPE (grid_real) :: meteoTemp !-------------------------end of declarations---------------------------------- !check if albedo and dem have been loaded by domain properties IF ( .NOT. dem_loaded) THEN CALL Catch ('error', 'SolarRadiation', & 'dem was not loaded ') END IF IF ( .NOT. albedo_loaded) THEN CALL Catch ('error', 'SolarRadiation', & 'albedo was not loaded ') END IF !check if air temperature and relative humidity have been initialized IF ( dtTemperature <= 0) THEN CALL Catch ('error', 'SolarRadiation', & 'air temperature not initialized ') END IF IF ( dtRelHumidity <= 0) THEN CALL Catch ('error', 'SolarRadiation', & 'air relative humidity not initialized ') END IF !set initial time timeNew = tstart !initialize grid CALL NewGrid (radiation, mask, 0.) CALL NewGrid (grid_devst, mask, 0.) !set deltat dtRadiation = IniReadInt ('dt', ini, section = 'solar-radiation') !check dt is multiple of dtMeteo IF (.NOT.(MOD(dtRadiation,dtMeteo) == 0)) THEN CALL Catch ('error', 'SolarRadiation', & 'dt is not multiple of dtMeteo') END IF !set valid threshold IF (KeyIsPresent ('valid-threshold', ini, & section = 'solar-radiation') ) THEN valid_prcn = IniReadReal ('valid-threshold', ini, & section = 'solar-radiation') ELSE ! set to default = 1.0 CALL Catch ('info', 'SolarRadiation', & 'valid-threshold not defined, set to 1') valid_prcn = 1. END IF !set cell-size cellsizeInterpolation = mask % cellsize !interpolation-assignment method IF (KeyIsPresent('interpolation-assignment', ini, & section = 'solar-radiation') ) THEN interpolationMethod_assignment = IniReadInt ('interpolation-assignment', & ini, section = 'solar-radiation') ELSE CALL Catch ('error', 'SolarRadiation', & 'interpolation-assignment missing in meteo configuration file') END IF !set interpolation method IF (interpolationMethod_assignment == 1) THEN interpolationMethod = IniReadInt ('interpolation', ini, & section = 'solar-radiation') CALL SetSpecificProperties (interpolationMethod, ini, mask) ELSE !read map interpolationMethod = -1 CALL GridByIni (ini, interpolationMethod_map, 'solar-radiation', & "interpolation") !scan for interpolation methods included interpolationMethod_vector = 0 DO i = 1, interpolationMethod_map % idim DO j = 1, interpolationMethod_map % jdim IF (interpolationMethod_map % mat(i,j) /= & interpolationMethod_map % nodata) THEN interpolationMethod_vector & (interpolationMethod_map % mat (i,j)) = & interpolationMethod_map % mat (i,j) END IF END DO END DO DO i = 1, 3 CALL SetSpecificProperties (interpolationMethod_vector (i), ini, mask) END DO END IF !allocate variable for netradiation CALL NewGrid (netRadiation, mask) CALL NewGrid (netRadiationFAO, mask) !scale factor and offset IF (KeyIsPresent('offset', ini, section = 'solar-radiation')) THEN offset_value = IniReadReal ('offset', ini, section = 'solar-radiation') ELSE offset_value = MISSING_DEF_REAL END IF IF (KeyIsPresent('scale-factor', ini, section = 'solar-radiation')) THEN scale_factor = IniReadReal ('scale-factor', ini, section = 'solar-radiation') ELSE scale_factor = MISSING_DEF_REAL END IF !set power_idw IF (KeyIsPresent('idw-power', ini, section = 'solar-radiation')) THEN idw_power = IniReadReal ('idw-power', ini, section = 'solar-radiation') ELSE !set default value idw_power = 2. END IF !file filename = IniReadString ('file', ini, section = 'solar-radiation') IF (interpolationMethod_assignment == 1 .AND. & interpolationMethod == 0) THEN !data are stored in net-cdf file !store net-cdf filename fileGrid = filename IF ( KeyIsPresent ('variable', ini, section = 'solar-radiation') ) THEN radiation % var_name = IniReadString ('variable', & ini, section = 'solar-radiation') !read grid and store as temporary file CALL NewGrid (layer = meteoTemp, fileName = TRIM(fileGrid), & fileFormat = NET_CDF, & variable = TRIM(radiation % var_name) ) ELSE IF (KeyIsPresent ('standard_name', ini, & section = 'solar-radiation') ) THEN radiation % standard_name = IniReadString ('standard_name', & ini, section = 'solar-radiation') ELSE CALL Catch ('error', 'SolarRadiation', & 'standard_name or variable missing in section temeprature' ) END IF !set cellsize to zero. Optimal cellsize is automatically computed cellsizeInterpolation = 0. ELSE !open file containing local measurements fileunit = GetUnit () OPEN(fileunit, file = filename(1:LEN_TRIM(filename)),status='old') END IF !use elevation to drift interpolation IF ( KeyIsPresent ('elevation-drift', ini, section = 'solar-radiation') ) THEN elevationDrift = IniReadInt ('elevation-drift', ini, section = & 'solar-radiation') ELSE elevationDrift = 0 !default, suppress drift END IF IF (elevationDrift == 1) THEN IF (interpolationMethod == 0 ) THEN CALL Catch ('error', 'SolarRadiation', & 'elevation drift cannot be applied when interpolation = 0 ') END IF IF ( KeyIsPresent ('time-zone', ini, section = 'solar-radiation') ) THEN timezone = IniReadInt ('time-zone', ini, section = & 'solar-radiation') ELSE CALL Catch ('error', 'SolarRadiation', & 'timezone is missing' ) END IF !set azimuth angles to compute shadow and convert to radians azimuth = (/0.0, 22.5, 45.0, 67.5, 90.0, 112.5, 135.0, 157.5, & 180.0, 202.5, 225.0, 247.5, 270.0, 292.5, 315.0, 337.5/) azimuth = azimuth * degToRad CALL DeriveSlope (dem, slope) CALL DeriveAspect (dem, aspect) CALL NewGrid (shadowGrid, dem) ALLOCATE ( viewangle (dem % idim,dem % jdim) ) END IF !grid exporting settings IF (KeyIsPresent ('export', ini, section = 'solar-radiation') ) THEN export = IniReadInt ('export', ini, section = 'solar-radiation') ELSE export = 0 END IF IF (export == 1) THEN !set export path name IF (KeyIsPresent ('export-path', ini, section = 'solar-radiation') ) THEN export_path = IniReadString ('export-path', ini, section = 'solar-radiation') !check if path ends with / or \ IF ( export_path (LEN_TRIM (export_path):LEN_TRIM (export_path)) /= '\' .AND. & export_path (LEN_TRIM (export_path):LEN_TRIM (export_path)) /= '/' ) THEN export_path (LEN_TRIM (export_path)+1:LEN_TRIM (export_path)+1) = '\' END IF IF (INDEX (export_path,'.') == 1) THEN !detected relative path, remove '.' export_path = export_path (2:LEN_TRIM(export_path)) !build absolute path export_path = TRIM(CurrentDir() ) // TRIM(export_path) END IF !check OS convention IF (GetOS () == WIN32) THEN export_path = ReplaceChar (export_path,'/','\') ELSE export_path = ReplaceChar (export_path,'\','/') END IF !check folder exists IF ( .NOT. DirExists (TRIM (export_path) ) ) THEN CALL Catch ('info', 'SolarRadiation', & 'creating directory: ', argument = TRIM(export_path)) CALL DirNew (export_path) END IF ELSE CALL Catch ('error', 'SolarRadiation', & 'missing export-path') END IF !starting time IF (KeyIsPresent ('export-start', ini, section = 'solar-radiation') ) THEN timeString = IniReadString ('export-start', ini, 'solar-radiation') export_start = timeString ELSE CALL Catch ('error', 'SolarRadiation', & 'missing export-start') END IF !initialize timeNewExport timeNewExport = export_start !ending time IF (KeyIsPresent ('export-stop', ini, section = 'solar-radiation') ) THEN timeString = IniReadString ('export-stop', ini, 'solar-radiation') export_stop = timeString ELSE CALL Catch ('error', 'SolarRadiation', & 'missing export-start') END IF !export dt IF (KeyIsPresent ('export-dt', ini, section = 'solar-radiation') ) THEN export_dt = IniReadInt ('export-dt', ini, section = 'solar-radiation') ELSE CALL Catch ('error', 'SolarRadiation', & 'missing export-dt') END IF !coordinate reference system IF (KeyIsPresent ('export-epsg', ini, section = 'solar-radiation') ) THEN export_epsg = IniReadInt ('export-epsg', ini, section = 'solar-radiation') exportGridMapping = DecodeEPSG (export_epsg) IF (exportGridMapping == radiation % grid_mapping) THEN needConversion = .FALSE. !initialize grid for exporting with CRS as meteo variable CALL NewGrid (exportedGrid, radiation) CALL NewGrid (exportedGridVar, radiation) CALL NewGrid (exportedGridNet, netRadiation) ELSE needConversion = .TRUE. !initialize grid for converting with required CRS exportedGrid % grid_mapping = DecodeEPSG (export_epsg) exportedGridVar % grid_mapping = DecodeEPSG (export_epsg) exportedGridNet % grid_mapping = DecodeEPSG (export_epsg) END IF ELSE CALL Catch ('info', 'SolarRadiation', & 'export-epsg not defined, use default') needConversion = .FALSE. !initialize grid for exporting with CRS CALL NewGrid (exportedGrid, radiation) CALL NewGrid (exportedGridVar, radiation) CALL NewGrid (exportedGridNet, netRadiation) END IF !export file format IF (KeyIsPresent ('export-format', ini, section = 'solar-radiation') ) THEN export_format = IniReadInt ('export-format', ini, section = 'solar-radiation') ELSE CALL Catch ('error', 'SolarRadiation', & 'missing export-format') END IF IF (export_format == 3) THEN !grid map CALL SetLongName ( 'downwelling_shortwave_flux_in_air', exportedGrid) CALL SetLongName ( 'net_downward_flux_in_air', exportedGridNet) CALL SetStandardName ( 'downwelling_shortwave_flux_in_air', exportedGrid) CALL SetStandardName ( 'net_downward_flux_in_air', exportedGridNet) CALL SetUnits ('Wm-2', exportedGrid) CALL SetUnits ('Wm-2', exportedGridNet) !if file exists, remove it export_file = TRIM(export_path) // 'radiation.nc' IF ( FileExists (export_file) ) THEN CALL FileDelete (export_file) END IF export_file_net = TRIM(export_path) // 'net_radiation.nc' IF ( FileExists (export_file_net) ) THEN CALL FileDelete (export_file_net) END IF !variance of kriging CALL SetLongName ( 'downwelling_shortwave_flux_in_air_variance', exportedGridVar) CALL SetStandardName ( 'downwelling_shortwave_flux_in_air_variance', exportedGridVar) CALL SetUnits ('Wm-2', exportedGridVar) !this unit is for exporting grid !if file exists, remove it export_file_var = TRIM(export_path) // 'radiation_variance.nc' IF ( FileExists (export_file_var) ) THEN CALL FileDelete (export_file_var) END IF END IF END IF !complete initialization IF (interpolationMethod == 0) THEN !Get the dt of imported field. Assume dt is regular dtGrid = GetDtGrid (filename = fileGrid, checkRegular = .TRUE.) !check dt is multiple of dtGrid IF (.NOT.(MOD(dtRadiation,dtGrid) == 0)) THEN CALL Catch ('error', 'SolarRadiation', & 'dt solar radiation is not multiple of dt of input grid') END IF ELSE !populate metadata CALL ReadMetadata (fileunit, radiometers) !check spatial reference system IF ( .NOT. radiometers % mapping == mask % grid_mapping) THEN CALL Catch ('info', 'SolarRadiation', & 'converting coordinate of stations') !convert stations' coordinate point1 % system = radiometers % mapping point2 % system = mask % grid_mapping radiometers % mapping = mask % grid_mapping DO i = 1, radiometers % countObs point1 % easting = radiometers % obs (i) % xyz % easting point1 % northing = radiometers % obs (i) % xyz % northing point1 % elevation = radiometers % obs (i) % z CALL Convert (point1, point2) radiometers % obs (i) % xyz % easting = point2 % easting radiometers % obs (i) % xyz % northing = point2 % northing END DO END IF END IF !interpolationMethod RETURN END SUBROUTINE SolarRadiationInit !============================================================================== !| Description: ! Read radiation data ! ! References: ! ! Ranzi R., Rosso R., Un modello idrologico distribuito, su base ! fisica, dello scioglimento nivale, Master thesis (in italian), ! Politecnico di Milano,1989. ! SUBROUTINE SolarRadiationRead & ! ( time, dem, albedo, temp, relHum ) IMPLICIT NONE !Arguments with intent(in): TYPE (DateTime), INTENT(IN) :: time !!current time TYPE (grid_real), INTENT(IN) :: dem !!used to apply drift of station data TYPE (grid_real), INTENT(IN) :: albedo !!used to apply drift of radiation site data TYPE (grid_real), INTENT(IN) :: temp !! air temperature (degree celsius) TYPE (grid_real), INTENT(IN) :: relHum !! air relative humidity (0-100) !local declarations: TYPE (DateTime) :: time_toread !! time to start reading from TYPE (DateTime) :: timeLocal !! local time for applying topographic drift CHARACTER (LEN = 300) :: filename CHARACTER (LEN = 300) :: varname REAL (KIND = float) :: rsquare INTEGER (KIND = short) :: i, j REAL (KIND = float) :: d !!solar declination REAL (KIND = float) :: w !!solar hour angle REAL (KIND = float) :: a1 !!sun elevation angle REAL (KIND = float) :: az !!azimuth REAL (KIND = float) :: mh !!atmosphere optical depth REAL (KIND = float) :: z !!terrain elevation REAL (KIND = float) :: Ic !!clear sky direct radiation REAL (KIND = float) :: D1 !!clear sky diffuse radiation REAL (KIND = float) :: DF !!diffuse radiation from surrounding elements REAL (KIND = float) :: A !!reflexed radiation from surrounding elements REAL (KIND = float) :: Rstar !! clear sky radiation = Ic + D1 REAL (KIND = float) :: radMeas !!measured value of radiation [W/m2] REAL (KIND = float) :: kt !! clearness index, cloudiness index complement REAL (KIND = float) :: fB1 !! hillslope coefficient REAL (KIND = float) :: cosT !! topographic factor REAL (KIND = float) :: Q !!direct radiation component on inclined surface REAL (KIND = float) :: radG !! ground radiation REAL (KIND = float), PARAMETER :: Isc = 1367. !! solar constant REAL (KIND = float), PARAMETER :: perclo = 0.22 !!minimum fraction of diffuse radiation REAL (KIND = float), PARAMETER :: k1 = 0.4 !!atmosphere diffusion coefficient !-------------------------end of declarations---------------------------------- IF ( .NOT. (time < timeNew) ) THEN !time_toread = time + - (dtRadiation - radiometers % timeIncrement) time_toread = time timeString = time_toread CALL Catch ('info', 'SolarRadiation', & 'read new solar radiation data: ', & argument = timeString) SELECT CASE (interpolationMethod) CASE (0) !input grid CALL ReadField (fileGrid, time_toread, & dtRadiation, dtGrid, & 'M', radiation, & varName = radiation % var_name) CASE DEFAULT !requires interpolation !read new data CALL ReadData (network = radiometers, fileunit = fileunit, & time = time_toread, aggr_time = dtRadiation, & aggr_type = 'ave', tresh = valid_prcn) IF (interpolationMethod_assignment == 1) THEN !only one method is applied CALL Interpolate (network = radiometers, & method = interpolationMethod, & near = neighbors, & idw_power = idw_power, & anisotropy = krige_anisotropy, & varmodel = krige_varmodel, & lags = krige_lags, & maxlag = krige_maxlag, & grid = radiation, & devst = grid_devst) ELSE !loop trough interpolation methods DO i = 1, 3 IF (interpolationMethod_vector(i) > 0) THEN CALL Interpolate (network = radiometers, & method = interpolationMethod_vector(i), & near = neighbors, & idw_power = idw_power, & anisotropy = krige_anisotropy, & varmodel = krige_varmodel, & lags = krige_lags, & maxlag = krige_maxlag, & grid = interpolatedMap (i), & devst = grid_devst) END IF END DO !compose the final interpolated field DO i = 1, interpolationMethod_map % idim DO j = 1, interpolationMethod_map % jdim IF (interpolationMethod_map % mat(i,j) /= & interpolationMethod_map % nodata ) THEN radiation % mat (i,j) = & interpolatedMap (interpolationMethod_map % mat(i,j)) % mat(i,j) END IF END DO END DO END IF !1 or more interpolation methods END SELECT IF (elevationDrift == 1) THEN !adjust downwelling shortwave !radiation to account for topography timeLocal = time + INT (timezone * hour) !compute solar declination d = SolarDeclination (timeLocal) !compute solar hour angle w = SolarHourAngle (timeLocal) !compute sun elevation angle a1 = SunElevationAngle (timeLocal, latCentroid) !compute azimuth az = SolarAzimuth (timeLocal, latCentroid) !compute shadow map CALL CastShadow (az, a1, viewangle, shadowGrid) DO i = 1, radiation % idim DO j = 1, radiation % jdim IF (radiation % mat (i,j) /= radiation % nodata ) THEN radMeas = radiation % mat (i,j) IF ( a1 <= 0 ) THEN !sun is below horizon radiation % mat (i,j) = radMeas netRadiation % mat (i,j) = radMeas netRadiationFAO % mat (i,j) = radMeas CYCLE END IF !elevation z = dem % mat (i,j) !optical depth mh = OpticalDepth (timeLocal, latCentroid, z) !clear sky direct radiation Ic = Isc * EXP ( -mh / SIN (a1) ) * SIN (a1) !clear sky diffuse radiation D1 = k1 * ( Isc * SIN (a1) - Ic ) !total clear sky radiation Rstar = Ic + D1 !clearness index, cloudiness factor complement (Ranzi, 1989) ! kt = 0 fully cloudy, kt = 1 clear sky kt = ( Rstar - radMeas ) / ( ( 1 - perclo ) * Rstar ) IF ( kt < 0. ) THEN kt = 0. END IF IF ( kt > 1. ) THEN kt = 1. END IF !diffuse radiation from surrounding elements IF ( Ic < 0. ) THEN DF = 0. ELSE DF = MIN (radMeas, D1 * ( 1 - Kt ) + radMeas * Kt) END IF !reflexed radiation from surrounding elements fB1 = 1. - slope % mat (i,j) A = radMeas * albedo % mat (i,j) * ( 1. - fB1 ) !direct radiation component cosT = COS (a1) * SIN (slope % mat(i,j) ) * & COS ( az - aspect % mat(i,j) ) + & SIN (a1) * COS (slope % mat(i,j) ) Q = ( radMeas - DF ) * cosT / SIN (a1) !check shadow IF ( shadowGrid % mat(i,j) == 1 .OR. Q < 0 ) THEN radG = DF + A cosT = -1 !error computing radG ELSE IF( ( Q + DF + A ) > Isc) THEN radG = Isc ELSE radG = Q + DF + A END IF radiation % mat (i,j) = radG !compute net radiation CALL ComputeNetRadiation ( kt, albedo % mat (i,j), & radiation % mat (i,j), temp % mat(i,j), & relHum % mat (i,j), netRadiation % mat (i,j) ) END IF END DO END DO ELSE !compute only cloudiness factor and net radiation timeLocal = time + INT (timezone * hour) !compute sun elevation angle a1 = SunElevationAngle (timeLocal, latCentroid) DO i = 1, radiation % idim DO j = 1, radiation % jdim IF (radiation % mat (i,j) /= radiation % nodata ) THEN radMeas = radiation % mat (i,j) IF ( a1 <= 0 ) THEN !sun is below horizon netRadiation % mat (i,j) = radMeas netRadiationFAO % mat (i,j) = radMeas CYCLE END IF !elevation z = dem % mat (i,j) !optical depth mh = OpticalDepth (timeLocal, latCentroid, z) !clear sky direct radiation Ic = Isc * EXP ( -mh / SIN (a1) ) * SIN (a1) !clear sky diffuse radiation D1 = k1 * ( Isc * SIN (a1) - Ic ) !total clear sky radiation Rstar = Ic + D1 !cloudiness factor complement. ! kt = 0 fully cloudy, kt = 1 clear sky kt = ( Rstar - radMeas ) / ( ( 1 - perclo ) * Rstar ) IF ( kt < 0. ) THEN kt = 0. END IF IF ( kt > 1. ) THEN kt = 1. END IF !compute net radiation CALL ComputeNetRadiation ( kt, albedo % mat (i,j), & radMeas, temp % mat(i,j), relHum % mat (i,j), & netRadiation % mat (i,j) ) !compute net radiation FAO with albedo = 0.23 CALL ComputeNetRadiation ( kt, 0.23, & radMeas, temp % mat(i,j), relHum % mat (i,j), & netRadiationFAO % mat (i,j) ) END IF END DO END DO END IF !elevationDrift ! END SELECT moved up !apply scale factor and offset, if defined IF (offset_value /= MISSING_DEF_REAL) THEN CALL Offset (radiation, offset_value) CALL Offset (netRadiation, offset_value) END IF IF (scale_factor /= MISSING_DEF_REAL) THEN CALL Scale (radiation, scale_factor) CALL Scale (netRadiation, scale_factor) END IF !grid exporting IF(export > 0 ) THEN IF( time == timeNewExport .AND. & time >= export_start .AND. & time <= export_stop ) THEN timeString = time timeString = timeString(1:19) timeString(14:14) = '-' timeString(17:17) = '-' !radiation grid_devst % reference_time = radiation % reference_time IF (needConversion) THEN CALL GridConvert (radiation, exportedGrid) CALL GridConvert (grid_devst, exportedGridVar) ELSE exportedGrid = radiation exportedGridVar = grid_devst END IF SELECT CASE (export_format) CASE (1) !esri-ascii export_file = TRIM(export_path) // TRIM(timeString) // & '_radiation.asc' CALL Catch ('info', 'SolarRadiation', & 'exporting grid to file: ', & argument = TRIM(export_file)) CALL ExportGrid (exportedGrid, export_file, ESRI_ASCII) IF (krige_var == 1) THEN !export kriging variance export_file_var = TRIM(export_path) // TRIM(timeString) // & '_radiation_variance.asc' CALL Catch ('info', 'SolarRadiation', & 'exporting variance grid to file: ', & argument = TRIM(export_file_var)) CALL ExportGrid (exportedGridVar, export_file_var, ESRI_ASCII) END IF CASE (2) !esri-binary export_file = TRIM(export_path) // TRIM(timeString) // & '_radiation' CALL Catch ('info', 'SolarRadiation', & 'exporting grid to file: ', & argument = TRIM(export_file)) CALL ExportGrid (exportedGrid, export_file, ESRI_BINARY) IF (krige_var == 1) THEN !export kriging variance export_file_var = TRIM(export_path) // TRIM(timeString) // & '_radiation_variance' CALL Catch ('info', 'SolarRadiation', & 'exporting variance grid to file: ', & argument = TRIM(export_file_var)) CALL ExportGrid (exportedGridVar, export_file_var, ESRI_BINARY) END IF CASE (3) !net_cdf CALL SetCurrentTime (time, exportedGrid) CALL Catch ('info', 'SolarRadiation', & 'exporting grid to file: ', & argument = TRIM(export_file)) CALL ExportGrid (exportedGrid, export_file, 'radiation', 'append') IF (krige_var == 1) THEN !export kriging variance CALL SetCurrentTime (time, exportedGridVar) CALL Catch ('info', 'SolarRadiation', & 'exporting grid to file: ', & argument = TRIM(export_file_var)) CALL ExportGrid (exportedGridVar, export_file_var, & 'radiation_variance', 'append') END IF END SELECT !net radiation IF (needConversion) THEN CALL GridConvert (netRadiation, exportedGridNet) ELSE exportedGridNet = netRadiation END IF SELECT CASE (export_format) CASE (1) !esri-ascii export_file = TRIM(export_path) // TRIM(timeString) // & '_net_radiation.asc' CALL Catch ('info', 'SolarRadiation', & 'exporting grid to file: ', & argument = TRIM(export_file)) CALL ExportGrid (exportedGridNet, export_file, ESRI_ASCII) CASE (2) !esri-binary export_file = TRIM(export_path) // TRIM(timeString) // & '_net_radiation' CALL Catch ('info', 'SolarRadiation', & 'exporting grid to file: ', & argument = TRIM(export_file)) CALL ExportGrid (exportedGridNet, export_file, ESRI_BINARY) CASE (3) !net_cdf CALL SetCurrentTime (time, exportedGridNet) CALL Catch ('info', 'SolarRadiation', & 'exporting grid to file: ', & argument = TRIM(export_file)) CALL ExportGrid (exportedGridNet, export_file_net, 'net_radiation', 'append') END SELECT timeNewExport = timeNewExport + export_dt END IF ENDIF !time forward timeNew = timeNew + dtRadiation END IF RETURN END SUBROUTINE SolarRadiationRead !============================================================================== !| Description: ! set properties and initialize variables for each interpolation method SUBROUTINE SetSpecificProperties & ! ( method, ini, mask ) IMPLICIT NONE !Arguments with intent(in): INTEGER (KIND = short), INTENT(in) :: method TYPE (IniList) , INTENT(in) :: ini TYPE (grid_integer) , INTENT(in) :: mask !----------------------------end of declarations------------------------------- !set data specific for each interpolation method SELECT CASE (method) CASE (1) !thiessen CALL NewGrid (interpolatedMap (1), mask, 0.) CASE (2) !inverse distance CALL NewGrid (interpolatedMap (2), mask, 0.) !nearest-points IF (KeyIsPresent ('nearest-points', ini, & section = 'solar-radiation') ) THEN neighbors = IniReadInt ('nearest-points', ini, & section = 'solar-radiation') ELSE !nearest-points missing CALL Catch ('error', 'SolarRadiation', & 'nearest-points missing in meteo configuration file') END IF CASE (3) !kriging CALL NewGrid (interpolatedMap (3), mask, 0.) !nearest-points IF (KeyIsPresent ('nearest-points', ini, & section = 'solar-radiation') ) THEN neighbors = IniReadInt ('nearest-points', ini, & section = 'solar-radiation') CALL Catch ('info', 'SolarRadiation', 'neighbors set to: ', & argument = ToString (neighbors) ) ELSE !nearest-points missing CALL Catch ('error', 'SolarRadiation', & 'nearest-points missing in meteo configuration file') END IF !kriging-variance IF (KeyIsPresent ('kriging-variance', ini, & section = 'solar-radiation') ) THEN krige_var = IniReadInt ('kriging-variance', ini, & section = 'solar-radiation') CALL Catch ('info', 'SolarRadiation', & 'kriging variance export set to ', & argument = ToString(krige_var) ) ELSE !kriging-variance missing krige_var = 0 CALL Catch ('info', 'SolarRadiation', & 'kriging variance export set to default (0)') END IF !kriging-anisotropy IF (KeyIsPresent ('kriging-anisotropy', ini, & section = 'solar-radiation') ) THEN krige_anisotropy = IniReadInt ('kriging-anisotropy', ini, & section = 'solar-radiation') CALL Catch ('info', 'SolarRadiation', 'kriging anisotropy set to: ', & argument = ToString (krige_anisotropy) ) ELSE !kriging-anisotropy missing krige_anisotropy = 0 CALL Catch ('info', 'SolarRadiation', & 'kriging anisotropy set to default (0)') END IF !kriging-variogram-model IF (KeyIsPresent ('kriging-variogram-model', ini, & section = 'solar-radiation') ) THEN krige_varmodel = IniReadInt ('kriging-variogram-model', ini, & section = 'solar-radiation') CALL Catch ('info', 'SolarRadiation', 'kriging variogram model set to: ', & argument = ToString (krige_varmodel) ) IF (krige_varmodel == 0 ) krige_varmodel = 5 !automatic fitting ELSE !kriging-variogram-model krige_varmodel = 2 CALL Catch ('info', 'SolarRadiation', & 'kriging variogram model set to default (2 exponential)') END IF !kriging-lags IF (KeyIsPresent ('kriging-lags', ini, & section = 'solar-radiation') ) THEN krige_lags = IniReadInt ('kriging-lags', ini, & section = 'solar-radiation') IF (krige_lags == 0) krige_lags = 15 CALL Catch ('info', 'SolarRadiation', 'kriging lags set to: ', & argument = ToString (krige_lags) ) ELSE !kriging-variogram-model krige_lags = 15 CALL Catch ('info', 'SolarRadiation', & 'kriging lags set to default (15)') END IF !kriging-maxlag IF (KeyIsPresent ('kriging-maxlag', ini, & section = 'solar-radiation') ) THEN krige_maxlag = IniReadInt ('kriging-maxlag', ini, & section = 'solar-radiation') CALL Catch ('info', 'SolarRadiation', 'kriging max lag set to: ', & argument = ToString (krige_maxlag) ) ELSE !kriging-maxlag missing krige_maxlag = 0 CALL Catch ('error', 'SolarRadiation', & 'kriging max lag set to default (0)') END IF END SELECT RETURN END SUBROUTINE SetSpecificProperties !============================================================================== !| Description: ! Compute the maximum angle of sky obstruction along 16 directions ! ! References: ! ! Zhang, Y., Chang, X. & Liang, J. Comparison of different algorithms for ! calculating the shading effects of topography on solar irradiance in ! a mountainous area. Environ Earth Sci 76, 295 (2017). ! https://doi.org/10.1007/s12665-017-6618-5 ! SUBROUTINE SkyView & ! ( dem, azimuth, view) IMPLICIT NONE !Arguments with intent(in): TYPE (grid_real) , INTENT(in) :: dem REAL (KIND = float), INTENT(in) :: azimuth (:) !azimuth angles [rad] !Arguments with intent(inout): TYPE (ViewingAngle), INTENT(inout) :: view(:,:) !local declarations INTEGER (KIND = short) :: i,j,l,i1,j1 TYPE (Coordinate) :: x0y0, x1y1 REAL (KIND = float) :: delta, deltax, deltay REAL (KIND = float) :: length LOGICAL :: check REAL (KIND = float) :: topAngle, topAngleMax !----------------------------end of declarations------------------------------- x0y0 % system = dem % grid_mapping x1y1 % system = dem % grid_mapping delta = dem % cellsize / 2. DO i = 1, dem % idim DO j = 1, dem % jdim IF (dem % mat (i,j) /= dem % nodata) THEN DO l = 1, SIZE (azimuth) !get coordinate of starting point CALL GetXY(i, j, dem, x0y0 % easting, x0y0 % northing) deltax = 0. deltay = 0. topAngle = 0. topAngleMax = 0. DO !add increment in x and y deltax = deltax + delta * SIN (azimuth (l) ) deltay = deltay + delta * COS (azimuth (l) ) x1y1 % easting = x0y0 % easting + deltax x1y1 % northing = x0y0 % northing + deltay !retrieve local coordinate of new cell CALL GetIJ (x1y1 % easting, x1y1 % northing, & dem, i1, j1, check) IF (.NOT. check ) THEN EXIT ELSE IF (dem % mat (i1, j1) == dem % nodata ) THEN EXIT END IF !compute distance length = Distance (x1y1, x0y0) !compute topographic angle topAngle = ATAN ( (dem % mat (i1, j1) - dem % mat (i,j) ) & / length ) IF (topAngle > topAngleMax) topAngleMax = topAngle END DO view(i,j) % angle (l) = topAngleMax END DO END IF END DO END DO RETURN END SUBROUTINE SkyView !============================================================================== !| Description: ! Compute shadow grid SUBROUTINE CastShadow & ! ( az, sunHeight, view, grid) IMPLICIT NONE !Arguments with intent(in): REAL (KIND = float), INTENT(in) :: az !current sun azimuth [rad] REAL (KIND = float), INTENT(in) :: sunHeight !current sun height [rad] TYPE (ViewingAngle), INTENT(in) :: view (:,:) !Arguments with intent(inout): TYPE (grid_real) , INTENT(inout) :: grid !local declarations INTEGER (KIND = short) :: i, j, m, l REAL (KIND = float) :: azDiff, azDiffMin !----------------------------end of declarations------------------------------- !search for closer azimuth to current sun azimuth azDiffMin = 100. l = 0 DO m = 1, SIZE (azimuth) azDiff = ABS ( azimuth (m) - az) IF ( azDiff < azDiffMin) THEN azDiffMin = azDiff l = m END IF END DO DO i = 1, grid % idim DO j = 1, grid % jdim IF ( grid % mat (i,j) /= grid % nodata ) THEN IF (sunHeight < view (i,j) % angle (l) ) THEN grid % mat (i,j) = 1 ELSE grid % mat (i,j) = 0 END IF END IF END DO END DO RETURN END SUBROUTINE CastShadow !============================================================================== !| Description: ! Compute solar declination. The declination of the sun is the angle ! between the equator and a line drawn from the centre of the Earth to ! the centre of the sun. ! ! References: ! ! Iqbal, M.: An introduction to solar radiation, Academis Press ! Canada, Ontario, 1983. FUNCTION SolarDeclination & ! (time) & ! RESULT (sdec) IMPLICIT NONE !Arguments with intent(in): TYPE (DateTime), INTENT(in) :: time !local declarations: REAL (KIND = float) :: sdec ![radians] REAL (KIND = float) :: dayAngle ![radians] INTEGER (KIND = short) :: dn !day number of the year !------------------------------------end of declarations----------------------- !day of year dn = DayOfYear (time, 'noleap') !day angle dayAngle = 2. * pi * (dn -1) / 365. !declination sdec = (0.006918 - 0.399912 * COS (dayAngle) + 0.070257 * SIN (dayAngle) - & 0.006758 * COS (2 * dayAngle) + 0.000907 * SIN (2 * dayAngle) - & 0.002697 * COS (3 * dayAngle) + 0.00148 * SIN (3 * dayAngle)) RETURN END FUNCTION SolarDeclination !============================================================================== !| Description: ! Compute the sun elevation angle ! ! References: ! ! Gates DM (1980) Biophysical ecology. Springer, New York, ! page 101, eq. 6.6 FUNCTION SunElevationAngle & ! (time, lat) & ! RESULT (sea) IMPLICIT NONE !Arguments with intent(in): TYPE (DateTime), INTENT(in) :: time REAL (KIND = float), INTENT(in) :: lat !latitude [radians] !local declarations: REAL (KIND = float) :: sea ! sun elevation angle [radians] REAL (KIND = float) :: dec !solar declination [radians] REAL (KIND = float) :: w !hour angle [radians] !------------------------------------end of declarations----------------------- w = SolarHourAngle (time) dec = SolarDeclination (time) sea = ASIN ( COS (lat) * COS (dec) * COS (w) + SIN (lat) * SIN (dec) ) RETURN END FUNCTION SunElevationAngle !============================================================================== !| Description: ! Compute the hour angle. the solar hour angle is an expression of time, ! hour angle is 0.000 degree, with the time before solar noon ! expressed as negative degrees, and the local time after solar noon ! expressed as positive degrees. For example, at 10:30 AM local ! apparent time the hour angle is -22.5 degree ! (15 degree per hour times 1.5 hours before noon). ! ! References: ! ! Iqbal, M.: An introduction to solar radiation, Academis Press ! Canada, Ontario, 1983. FUNCTION SolarHourAngle & ! (time) & ! RESULT (w) IMPLICIT NONE !Arguments with intent(in): TYPE (DateTime), INTENT(in) :: time !local declarations: REAL (KIND = float) :: w ! hour angle [radians] INTEGER (KIND = short) :: t !hour !------------------------------------end of declarations----------------------- t = GetHour (time) w = 15. * pi / 180. * (t - 12.) RETURN END FUNCTION SolarHourAngle !============================================================================== !| Description: ! Compute tthe atmospheric optical depth ! ! References: ! ! Kreith, F., and Kreider, J. F., 1978, Principles of Solar Engineering ! (New York: McGraw-Hill ). ! ! Cartwright, T. J., 1993, Modelling the World in a Spreadsheet: ! Environmental Simulation on a Microcomputer ! (Baltimore: Johns Hopkins University Press). ! ! LALIT KUMAR, ANDREW K. SKIDMORE & EDMUND KNOWLES (1997) Modelling ! topographic variation in solar radiation in a GIS environment, ! International Journal of Geographical Information Science, ! 11 :5, 475-497, DOI: 10.1080/136588197242266 FUNCTION OpticalDepth & ! (time, lat, z) & ! RESULT (s) IMPLICIT NONE !Arguments with intent(in): TYPE (DateTime), INTENT(in) :: time REAL (KIND = float), INTENT(in) :: lat !latitude [radians] REAL (KIND = float), INTENT(in) :: z !terrain elevation [m] !local declarations: REAL (KIND = float) :: s !!optical depth REAL (KIND = float) :: presscorr !atmospheric correction factor REAL (KIND = float) :: s0 !depth at sea level REAL (KIND = float) :: a1 ! sun elevation angle [radians] !------------------------------------end of declarations----------------------- !compute sun elevation angle a1 = SunElevationAngle (time, lat) !compute atmospheric optical depth at sea level s0 = ( 1229. + (614. * SIN (a1) ) **2. )**0.5 -614. * SIN (a1) !compute pressure correction factor presscorr = ( ( 288.15 - 0.0065 * z ) / 288.15 ) ** 5.25588 !optical depth s = s0 * presscorr RETURN END FUNCTION OpticalDepth !============================================================================== !| Description: ! Compute azimuth angle of the Sun's position in the north-clockwise ! convention ! ! References: ! ! Oke, T.R., Boundary layer climates, Second edition, Routledge, 1987. ! Appendix A1, eq. A1.2 FUNCTION SolarAzimuth & ! (time, lat) & ! RESULT (az) IMPLICIT NONE !Arguments with intent(in): TYPE (DateTime), INTENT(in) :: time REAL (KIND = float), INTENT(in) :: lat !latitude [radians] !local declarations: REAL (KIND = float) :: az !solar azimuth [radians] REAL (KIND = float) :: w ! hour angle [radians] REAL (KIND = float) :: d ! solar declination [radians] REAL (KIND = float) :: a1 ! sun elevation angle [radians] INTEGER (KIND = short) :: t !hour REAL (KIND = float) :: term !------------------------------------end of declarations----------------------- !compute hour angle w = SolarHourAngle (time) !compute declination d = SolarDeclination (time) !compute sun elevation angle a1 = SunElevationAngle (time, lat) !get hour t = GetHour (time) term = ( SIN (d) * COS (lat) - COS (d) * SIN (lat) * COS (w) ) / COS (a1) ! compute azimuth IF ( term < -1.) THEN az = pi ELSE IF (term > 1.) THEN az = 0. ELSE IF (t <= 12) THEN az = ACOS ( term ) ELSE az = 2. * pi - ACOS ( term ) END IF RETURN END FUNCTION SolarAzimuth !============================================================================== !| Description: ! Compute net radiation ! ! References: ! ! Wales-Smith, B. G. (1980). Estimates of net radiation for evaporation ! calculations, Hydrological Sciences Journal, 25:3, 237-242, ! DOI: 10.1080/02626668009491931 ! ! SUBROUTINE ComputeNetRadiation & ! (cfc, alb, short, hum, temp, netRad) IMPLICIT NONE !Arguments with intent(in): REAL (KIND = float), INTENT(in) :: cfc !!cloudiness factor complement REAL (KIND = float), INTENT(in) :: alb !!albedo REAL (KIND = float), INTENT(in) :: short !!shortwave radiation (W/m2) REAL (KIND = float), INTENT(in) :: hum !! air relative hunidity (0-100) REAL (KIND = float), INTENT(in) :: temp !! air temperature (degree celsius) !Arguments with intent(out): REAL (KIND = float), INTENT(out) :: netRad !!net radiation !local declarations: REAL (KIND = float) :: cf !!cloudiness factor REAL (KIND = float) :: netShort !!net shortwave radiation (W/m2) REAL (KIND = float) :: netLong !!net longwave radiation (W/m2) REAL (KIND = float) :: albedotemp !!temporary albedo value REAL (KIND = float) :: hum01 !!air relative humidity (0-1) REAL (KIND = float) :: humMin !!minimum value for air relative humidity (0-100) REAL (KIND = float) :: es !! saturation vapor pressure (Pa) REAL (KIND = float) :: ea !!actual vapor pressure (Pa) !------------------------------------end of declarations----------------------- !cloudiness factor cf = 1. - cfc !net shortwave radiation IF ( alb < 0. ) THEN albedotemp = 0.2 netShort = (1. - albedotemp) * short ELSE netShort = ( 1. - alb ) * short END IF !vapor pressure humMin = 10 IF ( hum > humMin ) THEN hum01 = hum / 100. ELSE hum01 = humMin / 100. END IF es = 6.107 * exp ( ( 17.27 * temp ) / ( temp + 237.3 ) ) ea = hum01 * es !net longwave radiation ( Wales-Smith, 1980) netLong = - ( stefanBoltzman * (temp + 273.15 )**4. ) * & ( 0.56 - 0.079 * ea**0.5 ) * ( 0.1 + 0.9 * cf ) !total net radiation netRad = netShort + netLong RETURN END SUBROUTINE ComputeNetRadiation END MODULE SolarRadiation